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Abstract 

In this work we use a low order dynamo model and study under which conditions can 
it reproduce solar grand minima. We begin by building the phase space of a proxy for the 
toroidal component of the solar magnetic field and we develop a model, derived from mean field 
dynamo theory, that gives the time evolution of the toroidal field. This model is characterized 
by a non-linear oscillator whose coefficients retain most of the physics behind dynamo theory. 
In the derivation of the model we also include stochastic oscillations in the a effect. We 
found no evidences that stochastic fluctuations in a linear a effect can trigger grand minima 
episodes in this model. In contrast, the model used points out that possible mechanism that 
can trigger grand minima should involve the meridional circulation, magnetic diffusivity or field 
intensification by buoyancy driven instabilities. 
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1. Introduction 

The Sun presents variability in several time scales, ranging from days to decades. The 
mechanisms behind this variability are still poorly understood although the common ground 
for most of them involve magnetic fields and turbulence. One of the main signatures of the solar 
magnetic activity is the cyclic formation of spots in the solar photosphere, usually known as 
sunspots. This sunspot cycle is also accompanied by changes in the solar spectrum. However, 
this cyclic activity is not regular since the peak amplitude and duration of the cycles changes 
with time. Sometimes these cycles even appear to be completely suppressed during long periods 
of time, giving rise to a specific kind of solar phenomena, the so called grand minima. In these 
periods the Sun appears to be in a very calm state, almost not exhibiting any sign of magnetic 
activity (spots, flares, etc.). The origin of these long periods of "solar inactivity" is still 
unknown and pose interesting scientific challenges. 

It is believed that the solar magnetic cycle has its origin in a dynamo process that operates 
in the convection zone and converts kinetic energy from the solar plasma flows into magnetic 
energy. When we have a grand minimum, the dynamo changes its operation regime and appar- 
ently shuts off for some time. The most famous grand minima that is registered is the Maunder 
Minima which occurred between the years of 1645 and 1715 (Eddy |8|). During this period, 
although there were no apparent signs of activity, several studies indicate that the dynamo was 
still operating (e.g. Beer, Tobias and Weiss [2|, Miyahara et al [10]). 
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To fully understand the intrinsic physics behind the dynamo one needs to resort to the 
magneto-hydrodynamic theory (MHD) which can be a very complex and difficult subject to 
fully grasp (Charbonneau 0]). Thankfully nowadays the fast development of computer science 
allows us to study these complex equations through the implementation of numerical dynamos. 
These "tools" , presently represent the best way of studying the processes involved in the dynamo 
operation. Some encouraging results on possible mechanisms behind grand minima have been 
presented in the last years (Charbonneau and Dikpati j^, Charbonneau, Blais-Laurier and 
St- Jean 15], Moss et al [12], Brandenburg and Spiegel [sl, Choudhuri and Karak [3]). 

As an alternative to MHD some authors, mainly in the 1990's, used low-dimensional chaotic 
systems to describe the behavior of the solar magnetic cycle (e.g. Ruzmaikin [isj, Ostriakov 
and Usoskin [l3|, Serre and Nesme-Ribes 1^). Low order models are simpler to compute but 
their interpretation can sometimes be tricky. Since they involve the collapse of the number of 
variables into a space with lower variables number, some information might be lost during the 
transformation. A low-order systems can be seen as a "projection" of a higher order system 
where the final result depends on the initial system and the "projection method" used. Due 
to this, it should be noted that reduced-order systems are often abstract representations which 
can loose physical meaning (Aantoulas and Sorensen [l|). By paying attention to these sensible 
points, dynamical system analysis involving low order models has proved to be a great tool in 



science. In more recent years, work developed by, e.g. Mininni et al Pontieri et al |l6lj . 
Wilmot-Smith i^l, Passos and Lopes suggests that within certain conditions, some of 

the observed properties of the solar magnetic field can be explained by low order dynamical 
models. 

In this work we intend to give a side perspective to the possible grand minima origins 
using a low-order dynamical system derived from dynamo theory. Although more limited than 
computational models this approach might be useful to build up intuition on physical processes. 

The model we use here is analogous to the one presented in Passos and Lopes [l3] and 
describes the evolution of the toroidal component of the solar magnetic field. Looking to the 
model's parameters we intend to study under which conditions can it reproduce grand minima. 
In order to compare this model with observational results, we use the sunspot number to build 
a proxy for the toroidal component and we look for the effects of grand minima in the phase 
space of this proxy. This gives us an experimental signature for grand minima that we should 
be able to reproduce with our model. We finish this present work with a discussion about the 
results obtained. 



2. Data and Grand Minima 

In order to study grand minima, we need to use solar activity records that go back in time 
to at least 1610, in order to include one of the most relevant grand minimum, the Maunder 
Minimum. For that purpose, we use the revised Sunspot Group Numbers (monthly averages), 
Rg, from Hoyt and Schatten \^ and available at NOAA database 0. After 1995 the time series 
is completed with the International Sunspot Number. 

As it is generally accepted, sunspots are a consequence of the toroidal magnetic field inside 
the convection zone, more specifically we can say that the sunspot number is proportional to 
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the magnetic energy (oc B"^) beneath the photosphere. Thus, we use Rg to build a proxy for this 
component of the field simply by assuming that B{t) oc ±\/Rg. To account for field reversals 
we change the sign of B{t) by hand for every sunspot cycle. To identify solar minima we used 
a low pass filter and selected the lowest values of the data series. Since identifying individual 
cycles in the period of the Maunder Minimum is very difficult, we decided to divide it into four 
separate "suppressed" cycles. At this point we would hke to note that since the amplitude of 
Rg during this period is very small, for the purpose of this work, a different choice would not 
have made an impact. In order to get the average behavior of the time series and eliminate 
"fast" transients (lower than 2.6 years), the proxy data is smoothed using a FFT filter (see 
figure (lH)). At this point we would like to note that the use of sunspots to build the B(t) proxy 
and the methodology applied, is going bind us to a characteristic dynamo scale whose behavior 
can, in principle, be reproduced by a low order model. 
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Figure 1: Top: Group Sunspot Number. Bottom: In black we have the built proxy for the toroidal field, B{t), 
superimposed to ±^/Rg in gray. 



As observed by Polygiannakis [I5l] a phase space reconstruction of the sunspot number 
hints that its behavior might be described by a non-linear oscillator. We pursue this idea but 
instead we use the proxy that we built. In order to construct our phase space, the numerical 
derivative, dB/dt, is computed using a time step of twelve months. 

Despite a small randomness, the trajectories of B{t) in the phase space appear to be stable, 
and seem to indicate that the solution for this oscillator is some kind of attractor. The only 
moment that the oscillator seems to seriously deviate from its "natural" action area (it collapses) 
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Figure 2: In this figure it's presented the constructed phase space of B{t) without the Maunder Minimum 
(right), and with it (left) in light gray dots to clearly distinguish it. 



is during the Maunder Minimum period (from approx. 1650 to 1720), depicted in gray in figure 
This is the experimental signature of grand minima that we will try to reproduce with the 
low-order model. 



3. Low order dynamo model with a stochastic ct effect 



In order to find an expression for a possible non-linear oscillator that might explain the 
behavior presented in the phase space of the toroidal field depicted in figure ([2]), we follow the 
ideas of Mininni et al llj] and Pontieri et al 16j. The model presented here is also discussed in 
Passos and Lopes ij] although with a different objective and derivation. Instead of proposing 
a purely mathematical inspired ad hoc expression for the oscillator, we intend to derive it from 
dynamo equations. This will allows us to connect physical mechanisms from the dynamo with 
coefficients in the oscillator's expression. 

We start by writing the equations for a mean field axisymmetric dynamo as shown in 
Charbonneau [6|. These equations give us the evolution of the mean solar magnetic field, B, 
classically decomposed into its toroidal and poloidal components, B = 
V X {ApC^). 



+ Bp with Bp 



dt 

OA, 
dt 



V 



+f [V X (ApC^)] ■ VVl 
V 



V2--U„-^-VffA 



(1) 

(2) 



where we have r = r sin 6, represents the differential rotation of the Sun, Vp is the flow in 
the meridional plane and t] is the magnetic diffusion. For simplification we will assume that r] 



is a constant in all of the convection zone { drj/dr = 0) and that the plasma is incompressible. 
We then get 



dt 



(3) 



dA 1 / 1 \ 

-gf = --:^rp-V{fAp) + {ao + ar{t))B^ + r]iv'--jAp, (4) 

where we introduced a simple linear a effect in the form of a = + Q;r(t) defined as having 
a constant part, Oq, and a stochastic part, ar-(^); that changes through time. The effect of 
a stochastic excitation in the a effect has been studied in numerical dynamo simulations and 
in theoretical dynamos by several authors (e.g. Charbonneau and Dikpati [3], Brandenburg 
and Spiegel [3:], Moss et al jl2|) and is justified by the angle dispersion around the mean tilt 
presented by bipolar active regions as they emerge. All these works present evidence that this 
stochastic effect might be behind grand minima phenomena, hence, we decided to introduce it. 
Also, in order to account for the removal of the toroidal field from the bottom of the convection 



zone by magnetic buoyancy we follow the suggestions of Pontieri et al [161] and add a term, 
r ~ jB'^/Snp, where 7 is a constant related to the buoyancy regime and p is the plasma 
density. This term can also work as an extra source for poloidal field due to buoyancy driven 



instabilities in the apex of rising flux tube (Rempel and Schussler [17|). As noted before, our 
sunspot derived proxy, gives us the magnetic field average behavior for a certain scale. In order 
to capture phenomena just on that scale, we truncate the dynamo equations by substituting 
V ^ I//0, where Iq is a specific length of interaction for the magnetic fields. 
After grouping terms in B^ and Ap we get 



dB^ 
dt 

dAr, 



dt 

where we have defined the coefficients, c„, as 



ciB^ + C2Ap - C3BI (5) 
CiAp + aoB^ + ar(t) B^, (6) 



^1 = (7) 



C2 = 7^ (8) 

We now concentrate in creating an expression for the time evolution of B^ since it is the 
quantity represented by our proxy B(t). To do so, we derive expression ([5]) in order to the time, 
and substitute (El) in it to take away the Ap dependence. After some mathematical manipulation 
it is possible to show that 

+ C2Mt))B, + p{3^Bl -l)^-XBl = 0, (10) 
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where uj"^ = — C2<yo, fJ^ = 2ci, ^ = C3/2C1 and A = C1C3 are coefficients that contain the solar 
physical structure (rotation, flows, diffusivity, etc.). 

The fact that the solar magnetic field presents a cyclic behavior in time, hints that the 
solution of this non-linear oscillator (van der Pol - Duffing type), in the phase space, would 
approximately correspond to a closed cycle (closed trajectory) or attractor whose shape depends 
on the structure parameters u, /i, ^ and A. From the dynamical system's point of view each one 
of these parameters will control the system in different ways. In the case without fluctuations in 
the a effect, i.e, ar{t) = 0, a; controls the frequency of the oscillations. In the phase space this 
coefficient is related to the amount of time the system takes to complete a closed trajectory, i.e., 
it can be defined asu = 2tt /T where T is the period of the cycle. The other term that also affects 
directly the frequency of oscillations is A. In the presence of a fluctuations, i.e. arif) 7^ 0, the 
term C2ar{t) will appear as a perturbation to the frequency. As for the remaining coefficients, 
/i controls the asymmetry between the rising and falling parts of the cycle and ^ affects directly 
the amplitude. 

In a ffist approximation, we can assume that the physical quantities involved in the structure 
parameters change in a time scale much longer than the magnetic field itself, e.g. from cycle 
to cycle, (like a quasi-static approach). Thus, if we look into what might change in these 
parameters that can create grand minima (a collapse to the center in terms of the phase space), 
one might derive some clues about the intervening physical mechanisms. 

3.1. Obtaining Minima with ar ? 

By analyzing equation (fTOj) we can see that the term (cj^ — C20ir{t)) is going to introduce 
variations in the system due to the random fluctuations of ar{t). For simplification lets call 
this term $r(^) and study its impact in the solution. 



$r.(t) = a;^ - C2ar-(t) = cj - C2ao - C2ar{t), 

|-^-^)'-?^(«o + Mt)), (11) 

^ V ' ' ' 

i) a) 

If one allows random ar{t) fluctuations to have values between ±ao then the solution's space 
for $r.(t) will contain an interval of values where the system is stable and another where the 
system is unstable. In the later interval the regular oscillations disappear and the solutions 
grow positive or negative depending on the value of the other coefficients and the point where 
the system is at that moment. Figure ([3]) shows a pictorial example of this. 



In this figure we used the average values for /i, ^ and A found in Passos and Lopes [14 
and present several scenarios for $r.(t). For the two black dashed trajectories, the system finds 
a stable solution after spending some time in an "almost saturated state". This happens for 
some super critical values of ^r{t)- On the other hand, for the same values but starting from a 
different point in the phase space (different initial conditions), the dark gray solid trajectories 
become unstable, drifting away from the attractor and growing indefinitely. A further example 
of the response of the system to a ^r(t) profile that, at a certain period, becomes super-critical 
and then recovers, is presented in figure (jl]). 
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Figure 3: Representation of the phase space (left) of equation (10) with /i = 0.1865, ^ = 0.0155, A = 0.0001 
and with ^r(t) — 0.1328 for four different initial conditions (dots). The light gray trajectory represents a stable 
solution with <f>r(i) — 0.2855 (corresponds to a period of 22 years). The time evolution of the B{t) for the same 
conditions is presented on the right. 



Even when we use a random ^r{t) profile we do not detect any collapse in the phase space, 
thus no minima. This can be caused by the fact that the buoyancy term, 7, might be working as 
a source keeping the dynamo alive when a weakens. At this point the only conclusion that we 
can derive is that, in this model, a linear a effect with random fluctuations cannot apparently 
be responsible for grand minima episodes. 

3.2. Obtaining Minima by maximizing ^ 

A way of making the trajectories in the phase space to collapse, is to increase the ^ parame- 
ter. Until a certain threshold, the higher ^ is, the lower is the amplitude of B{t). At this point, 
it is useful to remember that the construction of B{t) involves an unknown proportionality with 
Rg. This means that although in the phase space we built we can observe a full collapse of the 
trajectories, in order not to produce sunspots the field doesn't need to go all the way to zero. If 
the toroidal field strength falls bellow a certain threshold, then sunspot production will switch 
off. So, in principle, we could look into the physical processes contained in ,^ = c^b/icia + cib) 
and check which one can maximize this coefficient. 

e = . ^ (12) 

^ U '0 ; 

In this scenario three physical mechanisms come to play in the amplitude of the field: inten- 
sification of the field due to buoyancy instabilities (~ 7), magnetic diffusivity and meridional 
circulation amplitude. Several interconnections between these mechanisms can be arranged in 
order to increase ^. Since these mechanisms can also be affected by stochastic or other forcing 
factors (e.g. the meridional circulation can be affected by the field feedback, or 77 can incorpo- 
rate a quench) their inter-relations can create the necessary conditions to increase ^ and create 
a grand minimum. Figure ([5]) shows the response of the system to a temporal increase of ^. 
This increase makes B{t) collapse in the phase space, creating a signature analogous to a grand 
minima. 
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Figure 4: The evolution (top) and phase space (middle) of B{t) for a variable (bottom) that decreases 

to a super-critical value for some time and then recovers. The other coefficients remain constant and with the 
values used in figure 

4. Discussion 

Grand minima are episodes of solar "inactivity" that remain basically unexplained. Some 
clues about the physical mechanisms that might trigger these episodes of calmness in the Sun 
have been found through the use of computational solar dynamos. In this work we intend to 
present a different perspective on this subject by exploring a simple analytical model previously 
created to explain characteristics in the solar cycle. By analyzing the sunspot number time 
series, some information about the latest grand minima can be recovered. With this in mind, 
we use this time series to create an experimental proxy for the toroidal component of the solar 
magnetic field, B{t), construct its phase space and look for the signature of grand minima in 
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this phase space. The second step is to, based on dynamo theory, develop a simple model 
aimed to explain the average evolution to the toroidal field component. We modify the low- 



order model found in Passos and Lopes [ij] by adding stochastic fiuctuations to the a effect 



and we explore under which conditions can it reproduce physical solutions that resemble the 
grand minima signature found in the experimental data. 

As seen in section 3.1, in this simplified model, the presence of a stochastic component in a 
linear a effect can not produce a signature comparable to grand minima. Although this result 
is in apparent contrast with solutions found using more sophisticated computational dynamo 
models, the direct comparison is not trivial. In most of the computational models a non-linear 
a effect that incorporate quenching terms and/or different spacial locations is used. These 
non linearities introduce different physical dependencies in the structural coefficients, c„, and 
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ultimately a different final behavior of the system. We are currently working into incorporating 
such non-hnear effects in our model as well as other non-linearities in other coefficients. There 
is also the role of the buoyancy term, 7, that can act as source and whose interpretation is far 
from trivial. 

In this model, a "grand minima like" signature can be obtained by an increase of ^. In 
the example presented, the phase space docs not collapse all the way to the center, but this 
is no problem since studies of magnetic buoyancy of flux tubes have shown that there should 
be a threshold for in order for it to produce sunspots. The physical mechanisms that can 
contribute to the increase of $, are meridional velocity, magnetic diffusion and intensification 
of the field due to buoyancy instabilities. At this moment, trying to infer any other physical 
meaning from this result would be speculative. Nevertheless, it is our belief that the model can 
be improved in order to give more information about the interplay of dynamo mechanisms. For 
the time being, the model and method presented here should be seen as a proof of concept for 
an alternative way of tackling problems associated with the dynamo. 

We would like to thank the editor and the two anonymous referees for their constructing 
critics to the first version of this paper. They made us rethink about the problem and con- 
tributed to the overall improvement of the work. We also would like to thank the financial 
support of FCT (Fundagao para a Ciencia e Tecnologia, Portugal). 
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